Análisis Exploratorio de Eventos de Ocurrencia de Incendios

1. Introducción

Este notebook integra y compara la ocurrencia de incendios en la jurisdicción de la CAR a partir de tres fuentes:

• DGOAT CAR: Puntos críticos de incendios (Septiembre 2024).
• UNGRD: Emergencias registradas por la Unidad Nacional para la Gestión del Riesgo de Desastres.
• MODIS MCD64A1: Producto satelital de áreas quemadas.

Objetivos: 1. Cargar y limpiar cada fuente. 2. Calcular conteo de eventos y área anual por fuente. 3. Visualizar comparaciones temporales y espaciales.

2. Configuración del entorno y librerías

library(terra)
library(sf)
library(tidyverse)
library(readxl)
library(skimr)
library(leaflet)
library(RColorBrewer)
library(viridis)
library(lubridate)
library(patchwork)
library(plotly)

3. Comprobación de rutas de datos

# Rutas originales de archivos
data_paths <- c(
  "data/interm/base/jurisdiccion_car.gpkg",
  "data/interm/base/municipios.gpkg",
  "data/interm/eventos_ocurrencia/eventos_dgoat_geoambiental.gpkg",
  "data/interm/eventos_ocurrencia/eventos_dgoat_interno.csv",
  "data/interm/eventos_ocurrencia/eventos_ungrd_datosabiertos.csv",
  "data/interm/eventos_ocurrencia/eventos_ungrd_gestiondelriesgo.csv",
  "data/interm/satelital/modis_areas_quemadas.tif",
  "data/interm/eventos_ocurrencia/eventos_modis_mcd64a1.gpkg"
)
for (fp in data_paths) {
  if (!file.exists(fp)) stop(paste("Archivo no encontrado:", fp))
}

4. Carga de capas base

jurisdiccion_car <- terra::vect("data/interm/base/jurisdiccion_car.gpkg")
municipios_car    <- sf::read_sf("data/interm/base/municipios.gpkg")

plot(jurisdiccion_car, col = NA, lwd=2, border = "blue", main = "Jurisdicción CAR", )
plot(sf::st_geometry(municipios_car), add = TRUE, border = "darkgrey")

5. Carga de datos de eventos

# DGOAT
eventos_dgoat_geoambiental <- sf::read_sf("data/interm/eventos_ocurrencia/eventos_dgoat_geoambiental.gpkg")
eventos_dgoat_interno      <- readr::read_csv("data/interm/eventos_ocurrencia/eventos_dgoat_interno.csv")

# UNGRD
eventos_ungrd_datos_abiertos   <- readr::read_csv("data/interm/eventos_ocurrencia/eventos_ungrd_datosabiertos.csv")
eventos_ungrd_gestiondelriesgo <- readr::read_csv("data/interm/eventos_ocurrencia/eventos_ungrd_gestiondelriesgo.csv")

# MODIS
modis_areas_quemadas_raster <- terra::rast("data/interm/satelital/modis_areas_quemadas.tif")
eventos_modis_mcd64a1       <- sf::read_sf("data/interm/eventos_ocurrencia/eventos_modis_mcd64a1.gpkg")

6. Inspección rápida de datos

6.1 Registros DGOAT - Interno, Con Fecha Reporte

skimr::skim(eventos_dgoat_interno)

Data summary

Name	eventos_dgoat_interno
Number of rows	975
Number of columns	14
_______________________
Column type frequency:
character	7
Date	1
numeric	6
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
Municipio	0	1.00	4	26	0	101	0
Vereda Localidad	59	0.94	3	43	0	492	0
Nombre Punto	0	1.00	4	92	0	869	0
Clase	0	1.00	17	17	0	1	0
Causas Tipos	285	0.71	5	33	0	10	0
Categoria	250	0.74	12	40	0	8	0
Subcategoria	457	0.53	6	90	0	20	0

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
Fecha_reporte	177	0.82	2011-09-05	2023-11-01	2017-12-23	515

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Año	0	1.00	2017.82	2.87	2011.00	2016.00	2018.00	2019.00	2023	▂▅▇▆▃
Este	0	1.00	980069.53	41678.35	911183.00	946210.50	969155.00	1020314.50	1063265	▆▇▅▃▅
Norte	0	1.00	1032666.70	46354.45	916670.00	995299.00	1033648.00	1069699.00	1126979	▁▇▇▇▅
Altura	93	0.90	1700.85	945.84	2.71	915.00	1572.00	2613.50	3800	▅▆▃▇▁
area_metros	266	0.73	200372.31	522030.81	1.00	10400.00	37212.00	119900.00	5000000	▇▁▁▁▁
area_hectareas	266	0.73	20.04	52.20	0.00	1.04	3.72	11.99	500	▇▁▁▁▁

(eventos_dgoat_interno %>% group_by(Municipio) %>%
  summarise(n = n(), area_total_afectada = sum(area_hectareas, na.rm = T)) %>%
  arrange(desc(area_total_afectada)) %>% slice(1:20) %>%
  ggplot(aes(x = reorder(Municipio, area_total_afectada), y = area_total_afectada)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() + theme_minimal() +
  labs(title = "Area afectada por Municipio", x = "Municipio", y = "Area total afectada (ha)")) +

(eventos_dgoat_interno %>%
  count(Municipio) %>% arrange(desc(n)) %>% slice(1:20) %>%
  ggplot(aes(x = reorder(Municipio, n), y = n)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() + theme_minimal() +
  labs(title = "Registros por Municipio", x = "Municipio", y = "Cantidad de incendios"))

6.2 Registros UNGRD

skimr::skim(eventos_ungrd_gestiondelriesgo)

Data summary

Name	eventos_ungrd_gestiondelr…
Number of rows	3393
Number of columns	8
_______________________
Column type frequency:
character	5
Date	1
numeric	2
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
departamento	0	1	6	12	0	2	0
municipio	0	1	4	26	0	100	0
evento	0	1	8	29	0	4	0
observaciones	13	1	7	12461	0	3329	0
mes	0	1	3	4	0	12	0

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
fecha	0	1	2001-01-02	2024-12-16	2017-07-24	1372

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
area_hectareas	382	0.89	15.07	225.34	0	1	2	5	12000	▇▁▁▁▁
year	0	1.00	2017.79	4.68	2001	2015	2017	2023	2024	▁▁▃▇▆

(eventos_ungrd_gestiondelriesgo %>% group_by(municipio) %>%
  summarise(n = n(), area_total_afectada = sum(area_hectareas, na.rm = T)) %>%
  arrange(desc(area_total_afectada)) %>% slice(1:20) %>%
  ggplot(aes(x = reorder(municipio, area_total_afectada), y = area_total_afectada)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() + theme_minimal() +
  labs(title = "Area afectada por Municipio", x = "Municipio", y = "Area total afectada (ha)")) +

(eventos_ungrd_gestiondelriesgo %>%
  count(municipio) %>% arrange(desc(n)) %>% slice(1:20) %>%
  ggplot(aes(x = reorder(municipio, n), y = n)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() + theme_minimal() +
  labs(title = "Registros por Municipio", x = "Municipio", y = "Cantidad de incendios"))

6.3 Registros MODIS MCD64A1

skimr::skim(eventos_modis_mcd64a1)

Data summary

Name	eventos_modis_mcd64a1
Number of rows	145
Number of columns	7
_______________________
Column type frequency:
character	4
numeric	3
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
fecha	0	1	10	10	0	66	0
Municipio	0	1	4	19	0	39	0
CODDANE	0	1	5	5	0	39	0
geom	0	1	33	38	0	144	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
n_pixeles	0	1	7.55	9.23	1	2.00	4	9	55	▇▁▁▁▁
area_afectada	0	1	188.79	230.86	25	50.00	100	225	1375	▇▁▁▁▁
dia_juliano	0	1	194.70	90.64	11	165.67	228	247	363	▃▁▃▇▁

(eventos_modis_mcd64a1 %>% group_by(Municipio) %>%
  summarise(n = n(), area_total_afectada = sum(area_afectada, na.rm = T)) %>%
  arrange(desc(area_total_afectada)) %>% slice(1:20) %>%
  ggplot(aes(x = reorder(Municipio, area_total_afectada), y = area_total_afectada)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() + theme_minimal() +
  labs(title = "Area afectada por Municipio", x = "Municipio", y = "Area total afectada (ha)")) +

(eventos_modis_mcd64a1 %>%
  count(Municipio) %>% arrange(desc(n)) %>% slice(1:20) %>%
  ggplot(aes(x = reorder(Municipio, n), y = n)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() + theme_minimal() +
  labs(title = "Registros por Municipio", x = "Municipio", y = "Cantidad de incendios"))

7. Preprocesamiento y unificación

clean_municipio <- function(x) {
  stringi::stri_trans_general(x, "Latin-ASCII") %>% toupper()
}

# DGOAT
res_dgoat <- eventos_dgoat_interno %>%
  as_tibble() %>%
  mutate(
    Municipio = clean_municipio(Municipio),
    year      = as.numeric(Año)
  ) %>%
  filter(Clase == "Incendio Forestal") %>%
  group_by(Municipio, year) %>%
  summarise(
    n_dgoat    = n(),
    area_dgoat = sum(area_hectareas, na.rm = TRUE),
    .groups    = "drop"
  )

# UNGRD (portal gestión)
res_ungdr <- eventos_ungrd_gestiondelriesgo %>%
  as_tibble() %>%
  mutate(
    Municipio = clean_municipio(municipio),
    year      = year(fecha)
  ) %>%
  group_by(Municipio, year) %>%
  summarise(
    n_ungdr    = n(),
    area_ungdr = sum(area_hectareas, na.rm = TRUE),
    .groups    = "drop"
  )

# MODIS
res_modis <- eventos_modis_mcd64a1 %>%
  as_tibble() %>%
  mutate(
    Municipio = clean_municipio(Municipio),
    year      = year(fecha)
  ) %>%
  group_by(Municipio, year) %>%
  summarise(
    n_modis    = n(),
    area_modis = sum(area_afectada, na.rm = TRUE),
    .groups    = "drop"
  )

# Unión completa de fuentes
test_data_events <- res_dgoat %>%
  full_join(res_ungdr, by = c("Municipio", "year")) %>%
  full_join(res_modis, by = c("Municipio", "year"))

8. Resumen anual global

library(knitr)
resumen_por_year <- test_data_events %>%
  group_by(year) %>%
  summarise(
    total_n_dgoat    = sum(n_dgoat, na.rm = TRUE),
    total_n_ungdr    = sum(n_ungdr, na.rm = TRUE),
    total_n_modis    = sum(n_modis, na.rm = TRUE),
    total_area_dgoat = sum(area_dgoat, na.rm = TRUE),
    total_area_ungdr = sum(area_ungdr, na.rm = TRUE),
    total_area_modis = sum(area_modis, na.rm = TRUE),
    .groups = "drop"
  )
# Mostrar como tabla kable
kable(
  resumen_por_year,
  caption = "Resumen anual global: número de eventos y área quemada por fuente",
  col.names = c("Año", "Eventos DGOAT", "Eventos UNGRD", "Eventos MODIS", "Área DGOAT (ha)", "Área UNGRD (ha)", "Área MODIS (ha)"),
  format = "html",
  digits = 0
)

Resumen anual global: número de eventos y área quemada por fuente

Año	Eventos DGOAT	Eventos UNGRD	Eventos MODIS	Área DGOAT (ha)	Área UNGRD (ha)	Área MODIS (ha)
2001	0	54	0	0	476	0
2002	0	10	0	0	70	0
2003	0	15	0	0	224	0
2004	0	6	0	0	400	0
2005	0	1	0	0	600	0
2010	0	1	0	0	130	0
2011	4	0	10	5	0	2100
2012	30	167	12	337	2828	3400
2013	26	193	14	1028	1843	2725
2014	31	202	10	996	1843	2600
2015	142	512	22	3196	3909	6400
2016	122	458	4	2879	3016	1600
2017	91	194	1	590	515	25
2018	113	0	13	1386	0	1950
2019	177	405	6	1590	3758	975
2020	80	178	11	459	1720	1075
2021	29	74	5	77	329	400
2022	25	52	6	504	316	325
2023	105	287	9	1160	2956	1275
2024	0	584	22	0	20455	2525

9. Visualizaciones

9.1 Área quemada por fuente

df_area_long <- resumen_por_year %>%
  pivot_longer(
    cols        = starts_with("total_area"),
    names_to    = "Fuente",
    values_to   = "total_area",
    names_prefix= "total_area_"
  )

p_fuente <- ggplot(df_area_long, aes(x = year, y = total_area, color = Fuente)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::label_comma()) +
  labs(
    title    = "Área acumulada quemada por fuente",
    subtitle = "DGOAT vs UNGRD vs MODIS",
    x        = "Año",
    y        = "Área (ha)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(p_fuente)

9.2 Conteo de eventos anual

df_events_long <- resumen_por_year %>%
  pivot_longer(
    cols        = starts_with("total_n"),
    names_to    = "Fuente",
    values_to   = "total_events",
    names_prefix= "total_n_"
  )

p_count <- ggplot(df_events_long, aes(x = factor(year), y = total_events, fill = Fuente)) +
  geom_col(position = position_dodge()) +
  scale_y_continuous(labels = scales::label_comma()) +
  labs(
    title = "Conteo anual de eventos por fuente",
    x     = "Año",
    y     = "Número de eventos"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(p_count)

9.3 Mapa interactivo de incendios

10. Conclusiones y próximos pasos

• Se evidencian diferencias entre fuentes: DGOAT (detalle espacial), UNGRD (registro administrativo), MODIS (detección satelital).
  
• Validar superposición espacial y análisis de concordancia.
  
• Incorporar datos VIIRS para mayor resolución temporal.